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In this study the numerical performance of the fourth order compact formulation of the steady 
2-D incompressible Navier-Stokes equations introduced by Erturk et al. {Int. J. Numer. 
Methods Fluids, 50, 421-436) will be presented. The benchmark driven cavity flow problem 
will be solved using the introduced compact fourth order formulation of the Navier-Stokes 
equations with two different line iterative semi-implicit methods for both second and fourth 
order spatial accuracy. The extra CPU work needed for increasing the spatial accuracy from 

second order ((9( Ax^ )) to fourth order ((9( Ax"^ )) formulation will be presented. 
1. INTRODUCTION 

In Computational Fluid Dynamics (CFD) field of study High Order Compact 
formulations are becoming more popular. Compact formulations provide more 
accurate solutions in a compact stencil. 

In finite difference, in order to achieve fourth order spatial accuracy, standard five 
point discretization can be used. When a five point discretization is used, the points 
near the boundaries have to be treated specially. Another way to achieve fourth order 
spatial accuracy is to use High Order Compact schemes. High Order Compact 
schemes provide fourth order spatial accuracy in a 3x3 stencil, hence this type of 
formulations can be used near the boundaries without a complexity. 

In the literature, Zhang [14], Dennis and Hudson [3], MacKinnon and Johnson [6], 
Gupta et al. [5], Spotz and Carey [7] and Li et al. [9] have demonstrated the 
efficiency of high order compact schemes on the streamfunction and vorticity 
formulation of 2-D steady incompressible Navier-Stokes equations for uniform grids. 
Also in the literature the studies of Ge and Zhang [4] and Spotz and Carey [8] are 
example studies on the application of the high order compact scheme to nonuniform 
grids. The advantage of the high order compact schemes is that for a given flow 
problem and for a chosen grid mesh, high order compact formulation provides more 
accurate solutions (©(Ax"^)) compared to standard second order formulation 
((9( Ax^ )). Also for a given flow problem, the same level of accuracy of the solution 
obtained by standard second order formulation ((9(Ax^)) using a certain grid mesh, 
can be obtained with a smaller grid mesh when high order compact (©(Ax"^)) 
formulation is used. 

Recently Erturk and Gokcol [11] have presented a new fourth order compact 
formulation. The uniqueness of this formulation is that the presented "Compact 
Fourth Order Formulation of the Navier-Stokes" (EONS) equations are in the same 
form with the Navier-Stokes (NS) equations with additional coefficients. In fact the 
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NS equations are a subset of the PONS equations and obtained when the additional 
coefficients are chosen as zero. Therefore any numerical method that solves the NS 
equations can be easily applied to the introduced PONS equations to obtain fourth 
order spatial accurate solutions. Moreover, the most important feature of the PONS 
equations is that any existing code with second order accuracy ((9( Ax^ )) can easily be 
changed to provide fourth order accuracy ((9( Ax"^ )) by just adding some coefficients 
into the existing code at the expense of extra CPU work of evaluating these 
coefficients. This is an important feature such that if one already has a second order 
accurate code and wants to increase the accuracy of it to fourth order, instead of 
writing a new code, one can use the PONS equations and just by inserting some 
coefficient into the existing second order code, the existing second order code can 
turn into a fourth order code. Therefore the PONS equations introduced by Prturk and 
Gokcol [11] provide a very easy way to convert an existing second order accurate 
code into a fourth order accurate code and this is as simple as inserting some 
coefficients into the existing code. Of course, when this is done, the code will run a 
little slower because of the extra CPU work of evaluating the inserted coefficients. It 
will be good to estimate the CPU time needed for convergence of a converted fourth 
order accurate code compared to the CPU time needed for a second order accurate 
code. 

Zhang [15] have studied the convergence and performance of iterative methods with 
fourth order compact discretization schemes. To the best of the author's knowledge in 
the literature there is not a study that documents the numerical performance of high 
order compact formulation the Navier-Stokes equations compared to regular second 
order formulation of the Navier-Stokes equations in terms of numerical stability and 
convergence for a chosen iterative method. In this study using the PONS equations 
introduced by Erturk and Gokcol [11], we will numerically solve the Navier-Stokes 
equations for both fourth order (©(Ax"^)) and second order ((9(Ax^)) spatial 
accuracy. This way we will be able to compare the convergence and stability 
characteristics of both formulations. In this study we will also document the extra 
CPU work that is needed for convergence when a second order accurate code is 
converted into a fourth order accurate code using the introduced PONS equations by 
Erturk and Gokcol [11]. The stability and convergence characteristics of both 
formulations and also the extra CPU work can show variation depending on the 
iterative numerical method used for the solution therefore in this study we will use 
two different line iterative semi-implicit numerical methods. Using these two 
numerical methods we will solve the benchmark driven cavity flow problem. Pirst we 
will solve the cavity flow with second order ((9(Ax^)) spatial accuracy then we will 
solve the same flow with fourth order (OiAx^^)) spatial accuracy. We will document 
the stability characteristics, such as the maximum allowable time increment (AO, 
convergence characteristics, such as the number of iterations and the CPU time 
necessary for a chosen convergence criteria and also the extra CPU work that is 
needed to increase the spatial accuracy of the numerical solution from second order to 
fourth order using the PONS equations. 

2. FOURTH ORDER COMPACT FORMULATION 

In non-dimensional form, steady 2-D incompressible Navier-Stokes equations in 
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streamfunction {y/) and vorticity {co) formulation are given as 
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where x and y are the Cartesian coordinates and Re is the Reynolds number. 
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Erturk and Gokcol [11] have introduced the Fourth Order Navier-Stokes (PONS) 
equations. The introduced PONS equations are the following 
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As it is described briefly in Erturk and Gokcol [1 1], the numerical solutions of PONS 
equations (3) and (4) are fourth order accurate to NS equations (1) and (2), strictly 
provided that second order central discretizations shown in Table 1 are used and also 
strictly provided that a uniform grid mesh with Ax and Ay is used. We note that NS 
equations are a subset of PONS equations and obtained when the coefficients A, B , 
C , D , E and F in PONS equations are chosen as zero. The PONS equations are in 
the same form with the NS equations therefore any iterative numerical method applied 
to streamfunction and vorticity equation (1) and (2) can be easily applied to fourth 
order streamfunction and vorticity equation (3) and (4). Moreover if there is an 
existing code that solves the streamfunction and vorticity equation (1) and (2) with 
second order spatial accuracy, by just adding some coefficients A , B , C , D , E and 
F into the existing code using the PONS equations, the same existing code can 
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provide fourth order spatial accuracy. A single numerical code for the solution of the 
FONS equations can provide both second order and fourth order spatial accuracy by 
just setting some coefficients. Of course there is a pay off switching from second 
order to fourth order spatial accuracy, that is the extra cost of CPU work of 
calculating the coefficients A , B , C , D , E and F as defined in equation (5). 



3. FINITE DIFFERENCE EQUATIONS 



For numerical solutions of the Navier-Stokes equations (1) and (2), the following 
finite difference equations provide second order ((9( Ax:^ )) accuracy 
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where subscripts denote derivatives as defined in Table 1. 
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As explained in Erturk and Gokcol [11] the solution of the following finite difference 
equations are fourth order ((9( Ax:"^ )) accurate to the Navier-Stokes. 
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Note that equations (8) and (9) are in the same form with equations (6) and (7) except 
with additional coefficients A , B , C , D , E and F . In equations (8) and (9) if the 
coefficients A , B , C , D , E and F are chosen to be zero then the solution of these 
equations are second order accurate ((9( Ax', Ay')) to NS equations, since when the 
coefficients are zero, equations (8) and (9) are identical with equations (6) and (7). 
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However, in equations (8) and (9) if the coefficients A , B , C , D , E and F are 
calculated as they are defined in equation (10), then the solution of these equations 
(8) and (9) are fourth order accurate ((9( Ax"^, Ax^Ay^, Ay"^ )) to NS equations. When 
FONS equations are used, one can easily switch from second order to fourth order 
spatial accuracy using a single equation by just using the appropriate coefficients. 
Computationally, calculating the coefficients defined in equation (10) when fourth 
order accuracy is desired will require an extra CPU work compared to second order 
accuracy. In order to quantify the extra CPU work to switch from equations (6) and 
(7) to equations (8) and (9), i.e. switch from second order accuracy to fourth order 
accuracy, we solve the above equations with two different line iterative semi-implicit 
numerical method and document the CPU time for comparison. 

We note that the equations (8) and (9) are nonlinear equations therefore they need to 
be solved in an iterative manner. In order to have an iterative numerical algorithm we 
assign pseudo time derivatives to equations (8) and (9), thus we have 
dy/ 



■¥..+¥yy+0)-A (11) 



^ . ±(1 + 5)^^ + ±(1 + C)CD^ _ ^y, + D)CD^ + {y,^ + E)cD^ - F (12) 

ot Re Re 
We solve these equations (11) and (12) in the pseudo time domain until the solution 
converges to steady state. 

One of the numerical methods we will use to solve equations (11) and (12) is the 
Alternating Direction Implicit (ADI) method. ADI method is a very widely used 
numerical method and in this method a two dimensional problem is solved in two 
sweeps while solving the equation implicitly in one dimension in each sweep. The 
reader is referred to [12], [1] and [2] for details. When we apply the ADI method to 
solve equation (11), first we solve the following tri-diagonal system in x -direction 
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Similarly when we apply the ADI method to solve equation (12), we first solve the 
following tri-diagonal system in x -direction 
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where S^^ and S^^ denote the second order finite difference operators, and similarly 
and S denote the first order finite difference operators in x - and y -direction 



http://www.cavityflow.com 



respectively, for example 
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where / and j are the grid index and denote any differentiable quantity. 



The second numerical method we will use is the efficient numerical method proposed 
by Erturk et al. [10]. Following Erturk et al. [10], first using an implicit Euler 
approximation for the pseudo time derivatives in equations (11) and (12), we obtain 

the following finite difference formulations 
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We note that equations (18) and (19) are in fully implicit form and each equation 
requires the solution of a large banded matrix which is not computationally efficient. 
Instead, we spatially factorize these equations (18) and (19) thus we obtain the 
following finite difference equations 
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The advantage of these equations (20) and (21) is that each equation requires the 
solution of tridiagonal systems which is computationally very efficient using the 
Thomas algorithm. However, spatial factorization introduces l^f terms into the left 
hand side of equations (20) and (21), also these terms remain in the solution even at 
the steady state. To cancel out these l^f terms due to the factorization, Erturk et al, 
[10] have added the same amount of l^f terms to the right hand side of the equations 
so that the equations recover the correct physical representation at the steady state. 
The final form of the finite difference equations take the following form 
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M^{\ + C)5^^+M{y/^+E)5^co^ (23) 
y Ke J 

The reader is referred to Erturk et al. [10] for details. The solution methodology for 

equations (22) and (23) involves a two stage time level updating. For example, for the 

solution of equation (22), we first solve for the introduced variable / in x -direction 

in the following tri-diagonal system 
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When the solution for / is obtained, the streamfunction variable is advanced into the 
next time level by solving the following tri-diagonal system in y -direction 
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Similarly, for the solution of equation (23), we first solve for the introduced variable 
g in X -direction in the following tri-diagonal system 
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When the solution for g is obtained, the vorticity variable is advanced into the next 
time level by solving the following tri-diagonal system in y -direction 
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Stortkuhl et al. [13] have presented an analytical asymptotic solution near the comers 
of cavity and using finite element bilinear shape functions they also have presented a 
singularity removed boundary condition for vorticity at the comer points as well as at 
the wall points. For the boundary conditions, in both of the numerical methods 
described above we follow Stortkuhl et al. [13] and use the following expression for 
calculating vorticity values at the wall 
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where V is the speed of the wall which is equal to 1 for the moving top wall and 
equal to for the three stationary walls. 



For corner points, we again follow Stortkuhl et al. [13] and use the following 
expression for calculating the vorticity values 
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where again F is equal to 1 for tlie upper two comers and it is equal to for the 
bottom two comers. 



In explicit notation, for the wall points shown in Figure 1, the vorticity is calculated as 

the following 
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Similarly, for the comer points also shown in Figure 1, the vorticity is calculated as 

the following. 
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The reader is referred to Stortkuhl et al. [13] for details on the boundary conditions 



4. RESULTS 

In order to quantify the extra CPU work needed when a second order accuracy code is 
converted into a fourth order accuracy code using the FONS equations introduced by 
Erturk and Gokcol [11], we have solved both the second order and fourth order 
accurate equations (6), (7), (8) and (9) for the solution of the driven cavity flow. We 
consider the driven cavity flow for Reynolds numbers of i?e=100, 1000 and 3200, 
with using a grid mesh of 128 x 128 ( Ax = Ay = A/z ). We note that since we are mainly 
interested in finding the ratio of CPU time needed for convergence of a fourth order 
accuracy code to CPU time needed for convergence of a second order accuracy code, 
the choice of grid mesh size is not important, such that the ratio will be the same 
whether a coarse or fine grid mesh is used. 

In solving the equations we decided to use the Alternating Direction Implicit (ADI) 
method and the numerical method proposed by Erturk et al. [10]. By using two 
different numerical methods we would be able to see if the extra CPU work is 
dependent on the numerical method used. While doing this, as a by-product, we 
would also be able to compare the ADI method and the method proposed by Erturk et 
al. [10] in terms of numerical performance. In both of the numerical method we use, 
for both second order and fourth order accuracy, the two equations, i.e. the 
streamfunction and the vorticity equations, are solved separately. In order to 
document the extra CPU work when a fourth order accuracy is desired what we do is, 
we first solve for second order accuracy and solve equations (6) and (7). Then 
keeping the number of grids, the time step At and boundary conditions the same, we 
solve for fourth order accuracy thus we calculate and insert the coefficients A , B , C , 
D , E and F into the equations and solve for equations (8) and (9). While we solve 
the same flow problem, i.e. the driven cavity flow, for both second and fourth order 
accuracy we document the necessary number of iterations and the CPU time needed 
for a certain defined convergence criteria. This way we would be able to compare the 
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convergence characteristics of both second and fourth order formulations in terms of 
the number of iterations and the CPU time, and also we would be able to document 
the extra CPU time needed if a second order code is converted into a fourth order 
code using the PONS equations. 

For the choice of the time steps in solving the governing equations, we decided to use 
different time steps. A/ , for streamfunction and vorticity equations. In both of the 
numerical methods we use, while solving both the streamfunction and vorticity 
equations, tri-diagonal matrices appear on the implicit LHS of the equations. When 
second order accuracy is considered, in streamfunction equation the diagonal elements 
on the LHS matrices becomel + ^ , also in vorticity equations the diagonal elements 

on the LHS matrices become 1 + ^^. We choose different time steps for 

streamfunction and vorticity equations that would make the diagonal elements the 
same in both equations. Therefore for streamfunction equation we use = a/SJa^ and 
for vorticity equation we useA^ = aRel^h^ , where a is a coefficient we can choose. 
For fourth order accuracy we use the same time steps we use in second order 
accuracy. By using the same time steps in second and fourth order accuracy we would 
be able to compare the numerical stability of the PONS equations and the NS 
equations. In order to do this first we both solve the NS and the PONS equations 
using the same time steps. Then we increases, i.e. increase the time stepA^, and 
solve the NS and the PONS equations again. We continue doing this until at some A/ 
the solution does not converge. Therefore we would document the maximum 
allowable A/ for convergence for both the NS and the PONS equations for a given 
Reynolds number and grid mesh. This maximum allowable A^ for convergence is an 
indicative of the numerical stability. Por example, using either of the numerical 
method, i.e. the ADI method or the Erturk method, we solve the same flow problem 
using both second and fourth order formulations. Therefore for a chosen numerical 
method, the maximum allowable time step for second and fourth order formulations 
will be indicative of the numerical stability characteristics of the second order 
formulation compared to that of the fourth order formulation. Also, using either of the 
formulations, i.e. second or fourth order formulations, we solve the same flow 
problem using both the ADI method and the Erturk method. Therefore, for a chosen 
formulation, the maximum allowable time step for the ADI and the Erturk methods 
will be indicative of the numerical stability characteristics of the ADI method 
compared to that of the Erturk method. 

Our extensive numerical studies show that the increase in the extra CPU work is 
dependent on the computer and the compiler used. In this study, we run the codes on a 
64 bit HP ES45 machine with EV68 AlphaChip 1.25 GHz processors with HP Tru64 
UNIX operating system. We run the codes with both compiled normally and also 
compiled using maximum compiler optimization (-fast -05). 

We start the iterations from a homogenous initial guess and continue until a certain 
condition of convergence is satisfied. As the measure of the convergence to the same 
level, the residual of the equations can be used as it was also used in Erturk et al. [10]. 
However we are solving two different equations, the NS and the PONS equations, and 
try to compare the CPU time of convergence for each equation to the same level. 
Therefore the residual of these equations may not show the same convergence level. 
Alternatively, we can use the difference of the streamfunction and vorticity variables 
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between two time steps as the measure of convergence. However the solutions of the 
two different equations are shghtly different since one is spatially second order and 
the other is fourth order accurate. Since the solutions are different, the difference of 
the streamfunction and vorticity variables between two time steps may not also show 
the same convergence level for those equations. Therefore as the convergence criteria 
we decided to use the difference of the streamfunction and vorticity variables between 
two time steps normalized by the previous value of the corresponding variable, such 
that 
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These residuals provide an indication of the maximum percent change in y/ and co 
variables in each iteration step. In all of the data presented in this study, for both the 
solutions of the NS and the PONS equations, obtained using both of the numerical 
methods, we let the iterations converge until both Residual , and Residual , are less 

than 10"^ . At this convergence level, this would indicate that the variables y/ and co 
are changing less than 0.000001% of their value between two iterations at every grid 
point in the mesh. 



Figure 2, 3 and 4 show the streamline and vorticity contours of the driven cavity flow 
for 7?^ =100, 1000 and 3200 respectively, obtained using the method proposed by 
Erturk et al. [10] applied to PONS equations (©(Ax"^)). We note that both second 
order and fourth order accurate solutions of ADI method and the Erturk method, agree 
well with the solutions found in the literature especially with Erturk et al. [10], [11]. 

Using both of the numerical methods, we solve the driven cavity flow using different 
coefficients for time (a ), i.e. using different time steps, and document the CPU time 
and iteration number needed to the desired convergence level explained above. Table 
2 shows the CPU time and iteration numbers for ADI method for different Reynolds 
numbers using various a values. Table 3 shows the same for the method proposed by 
Erturk et al. [10]. Looking at Table 2 and 3, for both of the numerical methods the 
number of iterations for convergence is almost the same for second order and fourth 
order accuracy. However for both of the numerical methods the CPU time for fourth 
order accuracy is greater than the CPU time for second order accuracy as expected 
since the coefficients A , B , C , D , E and F have to be calculated at each iteration 
in fourth order accuracy which will result in an increase in the CPU time. The ratio of 
the CPU times of fourth order accuracy to second order accuracy show the increase in 
CPU time when we switch from second order accuracy to fourth order accuracy. From 
Tables 2 and 3 we see that this ratio seems to increase slightly when Reynolds number 
increases. 



It seems that for both of the numerical methods, at a given Reynolds number, the 
ratios of CPU time and iteration number for second and fourth order accuracy is 
constant and it is independent of the time step, i.e. a . 

For the ADI method, in Table 2, when the order of accuracy is increased to fourth 
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order from second order, the CPU time increases almost 1.76 times for 7?^ =100. This 
number increases as the Reynolds number increases and the increase in CPU time 
becomes 1.90 times for Re =3200. 

For the method proposed by Erturk et al. [10], in Table 3, when the order of accuracy 
is increased to fourth order from second order, the CPU time increases almost 1.60 
times for Re =100 and it is almost 1.68 times for Re =3200. 

We note that, when the order of accuracy is increased from second order to fourth 
order, the 1.6 and 1.68 times increase in CPU time for Re =100 and 3200 respectively 
for the method proposed by Erturk et al. [10] are less than the equivalent 1.76 and 
1.90 times increase in CPU time for the same Reynolds numbers for the ADI method. 
This shows that, the extra CPU time needed for fourth order accuracy when EONS 
equations are used, is dependent on the numerical method used and the extra CPU 
time for the method proposed by Erturk et al. [10] is much lower than the extra CPU 
time for the ADI method. 

In Table 2 and 3, comparing the two methods, for the same Reynolds numbers and for 
the same a values (AO, the iteration numbers for convergence are almost the same 
for both ADI and the method proposed by Erturk et al. [10], however the CPU time 
for ADI is less than that of the method proposed by Erturk et al. [10]. The reason for 
this is that in the method proposed by Erturk et al. [10] on the RHS of the finite 
difference equations more terms have to be calculated at each iteration and this 
increases the CPU time compared to the ADI method. 

For faster convergence one can use larger time steps, if the numerical method used 
has a higher numerical stability limit. Therefore, for a numerical method, the 
maximum allowable time step (AO for convergence gives an indication of the 
numerical stability limit of the numerical method. Since in this study we have used 
two different numerical methods for the same flow problem, we decided to compare 
the numerical stability limit of the two methods applied to both the NS and the FONS 
equations, by finding the maximum allowable time step for convergence for both 
numerical methods. In order to find the maximum allowable time step for 
convergence, for a given Reynolds number we solve the second and fourth order 
equations using both of the numerical methods several times while increasing a with 
0.01 increments each time, until the solution no longer converges. 

For the ADI method the maximum allowable a for convergence is 0.79 for second 
order accuracy and it is 0.78 for fourth order accuracy for Re=\QQQ. For the method 
proposed by Erturk et al. [10] the maximum a values was 1.89 for second order 
accuracy and it was 1.75 for fourth order accuracy. This would indicate that one can 
use much larger time steps in Erturk method compared to the ADI method, for 
example, for i?^=1000 the method proposed by Erturk et al. [10] allows to use 2.4 
times larger time step for the NS equations and 2.2 times larger time step for the 
FONS equations than the ADI method. From this we can conclude that the Erturk 
method has better numerical stability characteristics compared to the ADI method. 
When the maximum allowable time steps are used, the required CPU time for the 
method proposed by Erturk et al. [10] is almost 0.53 of the required CPU time for the 
ADI method. This means that the method proposed by Erturk et al. [10] converge 
almost twice faster than the ADI method when the maximum allowable time steps are 
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used. 



Comparing the numerical stability of the NS and the PONS equations, we see that for 
a chosen numerical method the PONS equations have slightly less stability limit than 
that of the NS equations. Por example, for the ADI method and for 7?e=1000 the 
value of 0.79 for maximum allowable a for convergence for second order accuracy 
drop down to 0.78 when fourth order accuracy is used. Also, for the Erturk method 
and for Re =1000 the maximum allowable a value of 1.89 for second order accuracy 
drop down to 1.75 if we switch to fourth order accuracy. This would indicate that, for 
fourth order formulations the maximum allowable time step for convergence is lower 
than the maximum allowable time step for convergence for second order 
formulations. 

We then decided to run the same codes compiled with using the maximum compiler 
optimization (-fast -05). Table 4 and 5 document the CPU time and iteration 
numbers when compiler optimization is used for the same conditions documented in 
Table 2 and 3 respectively. Comparing the numbers in Table 4 and 2 for ADI method, 
compiler optimization decreases the necessary CPU time for convergence about 0.71 
times for second order accuracy and about 0.51 times for fourth order accuracy. Also 
comparing the same in Table 5 and 3 for the method proposed by Erturk et al, [10], 
compiler optimization decreases the necessary CPU time for convergence about 0.54 
and 0.45 for second order and fourth order accuracy respectively. These numbers 
show that compiler optimization decreases the CPU time significantly and for the 
numerical method proposed by Erturk et al, [10] the codes almost run twice faster in 
terms of CPU time when compiler optimization is used. 

2. CONCLUSIONS 

The PONS equations introduced by Erturk and Gokcol [11] are in the same form of 
the NS equations, therefore any numerical method that solve the Navier-Stokes 
equations can be easily applied to the PONS equations in order to obtain fourth order 
accurate solutions ((9( Ax"^ )). One of the features of the introduced PONS equations is 
that any existing code that solve the Navier-Stokes equations with second order 
accuracy ((9(Ax:^)) can be easily altered to provide fourth order accuracy ((9(Ax:'^)) 
just by adding some coefficients into the existing code using the PONS equations. 
This way, the accuracy of any second order code can be easily increased to fourth 
order, however there is a pay off for this increased accuracy, that is the extra CPU 
time for calculating the added coefficients. 

In this study we have solved the NS equations and the PONS equations to document 
the extra CPU time necessary for convergence when an existing second order accurate 
code is altered to provide fourth order accurate solutions. Por this we have used two 
different numerical methods. We find that the extra CPU time is slightly dependent on 
the numerical method used. Por the ADI method to obtain fourth order accurate 
solutions of driven cavity flow, the CPU time increases 1.8 times compared to second 
order accurate solutions for 7?e=1000. Also for the numerical method proposed by 
Erturk et al. [10], with the cost of 1.64 times the CPU time necessary for second order 
accuracy, a fourth order accurate solutions can be obtained for i?^=1000, using the 
PONS equations. 
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The FONS equations introduced by Erturk and Gokcol [11] provide a very easy way 
of obtaining fourth order accurate solutions by just adding some coefficients into an 
existing second order accurate code, at the expense of a minor increase in the CPU 
time. 
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